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^ ■ Abstract 

\ Using Schrodinger Functional methods, we compute the non-perturbative renormal- 

' isation and renormahsation group running of several four-fermion operators, in the 

framework of lattice simulations with two dynamical Wilson quarks. Two classes 
^ \ of operators have been targeted: (i) those with left-left current structure and four 

?H ' propagating quark fields; (ii) all operators containing two static quarks. In both 

cases, only the parity-odd contributions have been considered, being the ones that 
renormalise multiplicatively. Our results, once combined with future simulations of 
the corresponding lattice hadronic matrix elements, may be used for the computa- 
tion of phenomenological quantities of interest, such as Bk and Bb (the latter also 
in the static limit). 



1 Introduction 



Hadronic matrix elements of four-fermion operators play an important role in the 
study of CP violation via CKM unitarity triangle analyses, as well as in the under- 
standing of the AI =1/2 enhancement puzzle in K ^ tttt decays. The only known 
technique to compute hadronic matrix elements from first principles, namely lat- 
tice QCD, has long been hampered by a number of systematic uncertainties. Most 
notably, the high computational cost of including light dynamical quarks in the 
simulations has enforced either the quenched approximation, or the use of heavy 
dynamical quark masses, which necessitate long and potentially uncontrolled ex- 
trapolations to the chiral regime. It is thus important to upgrade existing quenched 
results by the inclusion of dynamical fermion effects. For recent progress reports on 
lattice results on flavour Physics, see [1]. 

The present work is a step in this direction. The non-perturbative renormalisa- 
tion of four-fermion operators, as well as the corresponding renormalisation group 
(RG) running between hadronic scales of O(Aqcd) and perturbative ones of about 
100 GeV, is a necessary ingredient in the process of producing the properly renor- 
malised matrix element in the continuum limit. Quantities that determine these 
renormalisation properties have been computed in the quenched approximation, us- 
ing finite-size scaling techniques, for a broad class of four-fermion operators [2-4]. 
The regularisation of choice was that of non-perturbatively 0(a) improved Wilson 
quarks and standard plaquette gauge action; the renormalisation schemes used were 
of the Schrodinger functional (SF) type. The aim of the present work is to extend 
these results to QCD with N{ = 2 dynamical quarks. More specifically, we present 
results for: (i) the RG-running of left-left current relativistic four-fermion operators; 
(ii) the RG-running of all AB = 2 operators with two static heavy quarks; (iii) 
the renormalisation factors that match the above operators to their renormalisation 
group invariant (RGI) counterparts. The latter have been computed for a regular- 
isation of the relativistic quarks by the non-perturbatively 0(a) improved Wilson 
action. Preliminary results have been presented in [5]. 

As we will point out below, the results of this work are relevant for the com- 
putation of physical quantities such as the bag parameters Bk and Bb- In the 
quenched approximation, the combination of this renormalisation programme with 
computations of bare hadronic matrix elements has already produced high precision 
estimates of a few physical quantities in the continuum [6,7]. Moreover, knowledge 
of the continuum RGI operators, computed with Wilson fermions, has allowed the 
determination, through a matching procedure, of the renormalisation factors of the 
same operators in the Neuberger fermion regularisation [8]. 

The paper is organised as follows. In sect. 2 we introduce the operator basis 
and the renormalisation schemes adopted in the present work. We also recall some 
basic formulae used for the reconstruction of the operator scale evolution in the SF 
framework. Sect. 3 is devoted to a detailed description of the lattice simulations and 
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numerical analyses of the operator RG-running. In sect. 4 we discuss the renormal- 
isation of the four-quark operators at a low energy matching scale. Conclusions are 
drawn in sect. 5. In order to improve readability, some tables and figures have been 
collected at the end of the paper. 



2 Definitions and setup 

2.1 Renormalisation of four-fermion operators 

We will consider two different classes of four-fermion (dimension-six) operators: 

Or,r,ix) = \ [(V5i(x)riV'2(x)) (^3(a;)r2V'4(x)) ± (V'2 ^ ^^4)] , (2.1) 

^Sr^l^) = \ [(V^h(x)riV^i(x)) (4(x)r2V2(x)) ± (Vi ^^2)] . (2.2) 

In the above expressions Vfc is a relativistic quark field with flavour index /c, -i/'j^ 
(■0j^) are static (anti-)quark fields, F; are Dirac (spin) matrices, and the parentheses 
indicate summation over spin and colour indices. In the present formalism, all quark 
flavours are distinct, enabling us to separate the calculation of the scale-dependent 
logarithmic divergences, which is the aim of the present work, from the problem of 
eventual mixing with lower-dimensional operator^. 

The renormalisation pattern of the above operators is determined by the sym- 
metries of the regularised theory. In the the parity-odd sector, complete bases of 
operators in the relativistic and static cases are given by 

Qk ^ {^VA+AV ^VA-AV ^SP-PS' ^SP+PS' ^TT' } ' fc — 1,...,5 , (2-3) 

Q-k ^ {^vA+Av ^vA-AV) ^sp-ps! ^sp+ps' } ' A;=l,...,4, (2.4) 

respectively. The notation is standard and self-explanatory, indicating the operator 
spin matrices F/, with say, Ova+av = Ova + 0%^. A full analysis of the renor- 
malisation properties of these operator bases with relativistic Wilson fermions has 
been performed in [9, 10]. A result of these works which is of particular relevance is 
that, contrary to the parity-even case, characterised by operator mixing due to the 
explicit breaking of chiral symmetry by the Wilson term in the quark action, the 
parity-odd operators are protected by discrete symmetries, and hence their renor- 
malisation pattern is continuum-like [11]. We point out that RG-running is identical 
for parity-even and parity-odd operators of the same chiral representation, since in 
the continuum limit chiral symmetry is restored. On the other hand, the (physically 



^ These power subtractions typically appear for some specific choices of quark masses and/or 
flavour content (e.g. penguin operators). Their determination is independent of that of the loga- 
rithmic divergences, once mass independent renormalisation schemes are employed. 



2 



relevant) matrix elements of the parity-even operators can be mapped exactly to 
those of the parity-odd ones via the addition of a chirally twisted mass term to the 
lattice quark action [6,10,12]. 

From now on, we will consider the subset of operators 

Qt , Q'^g{Q^Q^ + 4Q+,Q+ + 2Q+,Q+-2Q+} . (2.5) 

All these operators renormalise multiplicatively; i.e. given an operator O e {Qf, Q'^} 
the corresponding operator insertion in any on-shell renormalised correlation func- 
tion is given by 

Oh,(x, /x) = lini Z{go,afi) 0{x; go) , (2.6) 

a— >0 

where go, a are the bare coupling and lattice spacing, respectively and fi is the renor- 
malisation scale. The RG-running of the operator is controlled by the anomalous 
dimension 7, defined by the Callan-Symanzik equation 

d 

IJ'-Q^On{x,fi) =7(5r(/i)) Oji{x,fi) , (2.7) 

supplemented by the corresponding RG-equation for the renormalised coupling g, 

/x|-g(^)=/3(g(/i)) . (2.8) 

In mass-independent renormalisation schemes, the /3-function and all anomalous 
dimensions depend only on the renormalised coupling g. They admit perturbative 
expansions of the form 

/9(5) (fco + + + • • •) , (2.9) 

7(5)''^° V(7o + 7i5' + 72/ + ...) , (2.10) 

in which the coefficients 60,61,70 are renormalisation scheme-independent. In par- 
ticular, the universal coefficients of the /?- function read 

and the universal leading order (LO) coefficients of the anomalous dimensions of the 
operators in Eq. <\2.5h are given by 

(2.12) 
(2.13) 
(2.14) 



70+ = 
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(47r)2 ' 
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Moreover, in the SF renormalisation scheme, the next-to-next-to-leading order (NNLO) 
coefficient of the /3-function is known to be [13] 



(47r)3 



0.483 - 0.275iVf + 0.0361iV/ - 0.00175iVf^ 



(2.15) 



Upon formal integration of Eq. (j2.7|) . one obtains the renormalisation group 
invariant (RGI) operator insertion 



OnGi{x) = 0^{x;ij) 



exp 



9 (m) 



dg 



f jial _ 7o_ 
\P{9) bog 



(2.16) 



while the RG evolution between two scales fJ.i,fJ-2 is given by the scaling factor 

Z{go,afi2) 



/_ 



lim 



) I3{g) j "-^0 Z{go,ani) ' 



(2.17) 



2.2 Schrodinger Functional renormalisation schemes 

Eqs. (I2.16p - (l2.17p are the starting point for the non-perturbative computation of 
the RG evolution of composite operators. To that purpose we introduce a family of 
Schrodinger Functional renormalisation schemes. The latter are defined by setting 
up the theory on a four-dimensional hypercube of physical size T x with Dirichlet 
boundary conditions in Euclidean time and periodic boundary conditions in the 
spatial directions, up to a phase 6. We refer the reader to [14] for an introduction 
to the SF setup. In the present work we always choose T = L and 6 = 0.5. We 
also assume that no background field is present. The renormalisation scale is set as 

^i = i/L. 

Renormalisation conditions are imposed on SF correlators, following [2-4]; for 
the sake of completeness, we briefly outline the method. We first introduce bilinear 
boundary sources projected to zero external momentum, 

S,,,jr]=a6^6,(x)rC.,(y) , (2.18) 
xy 

Here T denotes a Dirac matrix, the flavour indices si^2 can assume both relativistic 
and static values and the fields C (CO represent functional derivatives with respect 
to the fermionic boundary fields of the SF at the initial (final) time xq = {xq = T). 
The four-quark operators are then treated as local insertions in the bulk of the SF, 



4 



Figure 1: Diagrammatic representation of correlation functions: relativistic four-quark cor- 
relators (first diagram from left), static four-quark correlators JF*^*^-* (third diagram from 
left), relativistic boundary-to-boundary correlators fl'^, (second diagram from left) and 
static boundary-to-boundary correlators fl^^ and (fourth diagram from left). Euclidean 
time goes from left to right. Single (double) lines represent relativistic (static) valence 
quarks. 

giving rise to the correlation functions 

^[rA,rB,rc](^o) = ^"'(^aird Qf (^) S2i[rA] S45[rB]) , (2.20) 

4'l,rB,rc](^o) = ^"'<^3h[rc] Q!t{x) 5^ih[rA] S23[rB]) . (2.21) 

Clearly, the Dirac matrices of the boundary sources must be chosen so that the 
correlators do not vanish trivially (e.g. due to parity). 

In the above definitions a "spectator" light quark has been introduced with 
flavour s = 5 for and s = 3 for This quark field has no Wick contractions 

with the valence quarks of the operator insertion and propagates straight from the 
initial to the final time boundary. Its role is merely to allow for parity-even correla- 
tors of parity odd four-fermion insertions without the need of introducing non-zero 
external momenta. 

In order to isolate the operator ultraviolet divergences in Eqs. ()2.20p - (j2.2ip . 
one has to remove the boundary sources' additional divergences. To this end, we 
introduce a set of boundary-to-boundary correlators, 

/r^^ = -^(S;,.j75]S,,,j75]), (2.22) 

1 ^ 

= --^ Y.(KsM ^s,sAlk]) , (2.23) 
fc=i 

where the flavour indices si^2 may assume once again either relativistic or static val- 
ues. Wick contractions of four-quark and boundary-to-boundary correlation func- 
tions are depicted in Fig. [H 

Since the logarithmic divergences of the boundary fields are removed by multi- 
plicative renormalisation factors and Z^, it can be easily recognised that ratios 
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of correlators such as 



are free of boundary divergences for any choice of the parameter a. Thus, the 
operators of interest are renormaUsed through the conditions 

Z^{go,a/L)R^{T/2)=R^{T/2) , (2.26) 

30=0 

Z^^\go,a/L)Tl^^\T/2) =n^^\T/2) , (2.27) 

go=o 

where all correlation functions are to be evaluated at the chiral point. 

A crucial observation is that all renormalisation factors thus obtained are flavour- 
blind, in the sense that they remove the logarithmic divergences from any four- 
fermion operator of a given Dirac structure, irrespective of its specific flavour con- 
tent. For instance, in /S.H = 2 transitions, such as -B^(s) — -Bd(s) oscillations in the 
static limit, one identifies ipi = 4'2 = d{s) as a down (strange) quark, and tp'i = u 
as an up quark; for either flavour identification the same Z^^'^^s renormalise the 
corresponding operator. Similarly, in the relativistic quark case, the dimension-six 
operator, be it A5 = 2 (with ipi = ^jj^ = s and ■02 = "04 = oi' = 2 (with 
"01 = V'3 = ^ ^-iid il^2 = = d{s)), is renormalised by the same . Also in this case 
we note the presence of the spectator quark -05 = u in the renormalisation condi- 
tion. These renormalisation factors also remove the logarithmic divergences of other 
dimension-six operators with the same Dirac structure but different flavour content. 
For example, even if the renormalisation of some AF = 1 operators is only complete 
after power subtractions are taken into account, their logarithmic divergences are 
removed by the same and Z^'''^ factors. 

Another issue, related to the above discussion, is that here we are working with 
N{ = 2 dynamical flavours. In large- volume simulations of hadronic matrix elements, 
the latter would be naturally identified with the up and down sea quarks. However, 
most matrix elements of interest involve additional propagating physical flavour^. 
Thus, it becomes necessary to address the effects of partial quenching, in the context 
of operator renormalisation. 

Due to their flavour-blind nature, the renormalisation factors and RG running 
contained in this work account for the scale dependence of any matrix element of 
four-fermion operator, computed on an A^f = 2 sea, provided that partial quenching 



^An exception is that of AB = 2 oscillations of the Bd meson in the static limit. 
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does not generate extra scale-dependent mixing, which spoils multiplicative renor- 
malisation. This is indeed the case for AF = 2 meson oscillations, since in the chiral 
limit the relevant symmetries (CPS for Qf and CPT plus heavy quark spin symme- 
try for Q'^) remain valid and protect these operators from new counterterms. On 
the contrary, unphysical mixing, generated by partial quenching, may become an 
issue for AF = 1 hadronic decays, at least for penguin contributions [15]. We stress 
that this problem is not specific to the renormalisation schemes under consideration. 



2.3 Step-scaling functions 

The non-perturbative study of the RG-evolution of our composite operators is based 
on the step-scaling functions (SSFs) 



f \ r /r\ /r\ Z{go,a/{2L)) 

a[u) = hm.L,{u,a/L) , l^{u,a/L) — 



a^o Z{go,a/L) 



, (2.28) 



with Z e According to Eq. ([TT71) . the SSFs describe the operator Re- 

running between the scales /ii = 1/L and /i2 = 1/(2-^^). The choice of the ratio 
A^i//^2 = 2 as the smallest positive integer is made with the aim of minimising the 
effects of the ultraviolet cutoff at finite values of the latter. 

In practice, the SSFs are simulated at several values of the lattice spacing for 
fixed physical size (inverse renormalisation scale) L. The corresponding values of the 
inverse bare coupling j3 = Q/qq are indeed tuned by requiring that the renormalised 
SF coupling 5gp, and hence L, are kept constant. The critical mass mcr is obtained 
from the requirement that the PC AC mass vanishes in the 0(a) improved theory. 

Once (j{u) is computed at several different values of the squared gauge cou- 
pling u, it is possible to reconstruct the RG evolution factor C/(/Xpt, /ihad) between 
a hadronic scale //had) in the range of a few hundred MeV, and a perturbative one 
//pt, in the high-energy regime. This in turn leads to the computation of the RGI 
operator of Eq. (|2.16p . with controlled systematic uncertainties, by splitting the 
exponential on the rhs of Eq. ()2.16p . evaluated at /x = /Xhadi as follows: 



c(/i; 



had J 



"5^(/"had)' 


70 
2f>0 




47r 


exp < 


L 



= c(/ipt) C/(/Upt,/ihad) ■ (2.29) 

The second factor on the rhs is known non-perturbatively as a product of continuum 
SSFs cr(u); cf. Eq. p.l7p . The first factor can be safely computed at next-to- leading 
order (NLO) in perturbation theory, provided the scale /Xpt is high enough to render 
NNLO effects negligible. The full procedure for the construction of U (/^pt, /^had) has 
been introduced in [16] for the running quark mass in the quenched approximation. 
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a 


Qt 


75 


75 


75 





Qi 


7fc 


75 


7fc 


1 


Q't 


75 


75 


75 


1/2 


Q't 


75 


75 


75 





Q't 


75 


75 


75 





Q't 


75 


75 


75 






Table 1: Optimal renormalisation schemes for the various four-quark operators. 

and subsequently applied in several contexts (for a recent review, see [17]). The 
reader is referred to the above-mentioned works for a detailed description of the 
method. More specifically, since the present work concerns two-flavour QCD, we 
follow closely the work of ref. [18] on the running quark mass with A'^f = 2 dynamical 
quarks. 

Concerning the choice of [rA,rB,rc] and a, we observe that in our quenched 
studies [2,4] we have considered five possible non-trivial Dirac structures that pre- 
serve cubic symmetry at vanishing external momenta, i.e. 

[rA,rB,rc] = { [75,75,75], (^ijkhi,ij,ik], [75,7^,7^], 

[7fc,75,7fc], [7fc,7fe,75] } (2.30) 

(where a sum over repeated indices is understood), and various possible values of the 
a parameter, namely a = {0, 1/2} for the static operators and a = {0, 1, 3/2} for the 
relativistic ones. Not all of the resulting renormalisation schemes were equally well 
suited to our purposes: some of them were characterised by a RG running with a 
slow perturbative convergence at NLO; cf. refs. [3,10]. This rendered the matching 
of perturbative and non-perturbative running at Upt (cf. Eq. (j2.29p ) less reliable 
and the systematics hard to control; see sect. 3.3 below for more details. The same 
considerations are valid in the present case of two dynamical quarks. For the sake 
of brevity we will concentrate only in those schemes which have been found to be 
best behaved in the present unquenched study. These optimal schemes are specified 
in Table [TJ A complete account of our results in all schemes considered is available 
upon request. 

The non-universal two-loop coefficients of the anomalous dimensions Eq. (|2.10|) 
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for the operators of interest in the aforementioned optimal schemes read [3, 4] 



7+ = — — r 0.0828(48) + 0.03200(28)Arfl , (2.31) 

(47r)^ 

7" = _^_r_o.6880(24) + 0.12648(16)iVfl , (2.32) 
(47r)^ 

= 1-345(2) + 0.0008(2)iVf] , (2.33) 

7?^ = -^[-1-251(1) + 0.11637(8)Arf] , (2.34) 

7f^ = --^[-0.327(3) + 0.1211(2)Arf] , (2.35) 

7^^ = -7^2 [-0.146(1) + 0.06784(8)iVf] . (2.36) 



3 Non-perturbative computation of the RG running 

3.1 Simulations details 

Our simulations are based on the regularisation of relativistic quarks by the non- 
perturbatively 0(a) improved Wilson action, with the Sheikoleslami-Wohlert (SW) 
coefficient Cgw determined in [19]. Static quarks have been instead discretised as 
proposed in [20]. In particular, all results reported in this work refer to the so- 
called HYP2 action, i.e. the lattice static action of [21], with the standard parallel 
transporter J7(0, x) replaced by the temporal hypercubic link introduced in [22], 
and a choice of the smearing parameters (01,02,03) = (1.0,1.0,0.5). The latter 
minimises the quenched static self-energy, providing the largest exponential increase 
of the signal-to-noisc ratio in the static-quark propagator, when compared to the 
original Eichtcn-Hill action. The minimum of the static self-energy is shifted by 
internal quark-loops only at NLO in perturbation theory: such shift is thus expected 
to be relatively small. 

With the above prescriptions, the SSFs have been computed at six different 
values of the SF renormalised coupling, corresponding to six different physical lattice 
lengths L. For each physical volume three different values of the lattice spacing have 
been simulated, corresponding to lattices with L/a = 6, 8, 12 (and 2L/a = 12, 16, 24 
respectively) for the computation of Z{gQ,a/ L) (and Z {go,a/{2L))). 

The gauge configuration ensemble used in the present work (generated with 
Nf = 2 dynamical fermions) and the tuning of the lattice parameters (/3, k) have 
been taken over from [18]. All technical details concerning these dynamical fermion 
simulations are discussed in that work. The one technical aspect that makes a 
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significant difference in our case concerns the perturbative value of the boundary 
improvement coefficient Ct [23]. As pointed out in [18], the gauge configurations at 
the three weakest couplings have been produced using the one-loop perturbative 
estimate of Ct [23], except for (L/o = 6, /? = 7.5420) and (L/a = 8, /? = 7.7206). 
For these two cases and for the three stronger couplings, the two-loop value of Ct [24] 
has been used. We have enforced the same Ct values in the valence propagators. 
Comparison of the results of two different simulations, namely 

= 1.5031(25) , L/a = 6, /3 = 7.5000 , k = 0.1338150 , ct = one - loop , 
= 1.5078(44) , L/a = 6, /3 = 7.5420 , k = 0.1337050 , ct = two - loop . 

shows that the renormalisation factor Zp of the pseudoscalar density, analysed in 
[18], is subject to a relative 4 per mille variation, corresponding to a mild discrepancy 
of about 1.5cr with regards to our statistical uncertainty. Unfortunately, for the four- 
fermion correlation functions this remains true only for the operator while the 
Q\ and Q!\ cases show relative variations of the order of 1-2%, corresponding to 
differences of about 2.5cj with regards to the statistical precision. We expect that, for 
a given renormalised coupling ^gp, this effect diminishes at finer lattice resolutions 
ajL (i.e. closer to the continuum), while it becomes more pronounced in the strong 
coupling region, at constant L/a. In principle, this problem can be removed by 
performing all simulations with a two-loop estimate of Ct and/or a smaller resolution 
ajL. As further dynamical simulations are beyond the scope of the present work, we 
limit ourselves in stating that our results are subject to this ill-controlled systematic 
uncertainty, which, in view of the fact that the one-loop value of Ct is only used 
in the weak coupling region, is however not expected to be significant for our final 
results. In this respect, we have checked that the (final) overall result is unaffected 
when either ct is employed at this coupling. It is also encouraging that, as we 
will see below, including or discarding the L/a = 6 data-points in the continuum 
extrapolations does not alter the final results significantly. 

Numerical results are collected in Tables ISHlOl Statistical errors were computed 
by a jackknife analysis. The estimates of the autocorrelation times, calculated with 
the autocorrelation function method, the method of ref. [25] and the binning method, 
were found to be compatible. 

3.2 Continuum extrapolation of the step-scaling functions 

Since we do not implement 0(a) improvement of four-fermion operators, the only 
linear cutoff effects that are removed from S(u, a/L) are those cancelled by the SW 
term in the fermion action. Therefore, we expect SSFs to approach the continuum 
limit linearly in a/L and correspondingly we fit to the ansatz 
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0.9793 


1.010(11) 


0.983(07) 


0.946(08) 


1.004(07) 


0.960(05) 


0.990(05) 


1.1814 


1.044(15) 


0.965(10) 


0.951(12) 


0.991(08) 


0.942(07) 


0.976(05) 


1.5078 


1.039(21) 


0.953(11) 


0.932(13) 


0.987(10) 


0.932(09) 


0.970(07) 


2.0142 


1.040(18) 


0.936(11) 


0.896(11) 


0.985(10) 


0.901(09) 


0.955(08) 


2.4792 


1.078(35) 


0.879(19) 


0.890(19) 


0.958(13) 


0.873(14) 


0.938(12) 


3.3340 


1.129(37) 


0.862(25) 


0.784(23) 


0.938(17) 


0.798(18) 


0.905(16) 



Table 2: Results of the continuum limit extrapolation of the lattice step-scaling functions 
and S^*^^ . Data have been fitted from all available lattice resolutions as linear functions 
in [a/L). 

In practice it is often observed that the data corresponding to L/a = 8, 12 are 
compatible within errors, whereas the L/a = Q result, bearing the largest cutoff 
effects, is off. This suggests that, in analogy to [18], a weighted average of the two 
finest lattice results may be a reliable estimate of the continuum limit value. We 
have checked that, in most cases, linear fits to all three data-points and weighted 
averages of the two results from the finer lattices lead to continuum limit estimates, 
compatible within one standard deviation; cf. Figs. [21 HI and El Fit results are 
reported in Table [2l Since the discretisation errors are 0(a) and not O(a^) as 
in [18], we conservatively quote, as our best results, those obtained from linear 
extrapolations involving all three data-points. 

It should be added that, besides the HYP2 action, we have also tried other static 
quark action varieties, namely the APE and the HYPl ones (see ref. [20]), which 
differ from HYP2 by 0{a?) lattice artefacts. Since the static four-quark operators 
are not 0(o) improved, it is reasonable to expect significant discretisation effects at 
the coarsest lattice spacing, which would enable combined fits of the data from all 
three actions, constrained to a common value in the continuum limit. Unsurpris- 
ingly, the situation turned out to be similar to that of [4] , in that data obtained with 
the above actions do not differ noticeably (even at L/a = 6) and are very strongly 
correlated. Consequently, a combined continuum extrapolation affects the contin- 
uum limit results only marginally, with the relative error decreasing only by a few 
percent. For this reason we only quote results from the HYP 2 analysis. 

3.3 RG running in the continuum limit 

In order to compute the RG running of the operators in the continuum limit as 
described in [18], the continuum SSFs have to be fitted to some functional form. 
The simplest choice is represented by a polynomial 

a{u) = l + siu + S2U^ + s^u^ + ... , (3.2) 
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n 


u 


\ in ax / 




c(^'(L-i ) 


c(2'(L-i ) 


s(3)|'r-l \ 

\ max / 


A(4)|'r-l ^ 

v niax / 





4.610 


1.246 


0.551 


0.807 


0.680 


0.524 


0.776 


1 


3.032 


1.225(26) 


0.554(10) 


0.775(14) 


0.732(09) 


0.532(08) 


0.788(09) 


o 
2 


z.o4i 


1 oi o/oo^ 
L.212[6o) 


0.5DD(14) 


0. / /3(2U) 


0. (51(12 ) 


0.538(11 ) 


0. (94(12) 


3 


1.918 


1.205(46) 


0.564(16) 


0.773(24) 


0.759(15) 


0.541(13) 


0.797(14) 


4 


1.628 


1.202(53) 


0.561(18) 


0.772(27) 


0.762(18) 


0.541(14) 


0.797(16) 


5 


1.414 


1.201(60) 


0.558(20) 


0.772(30) 


0.763(20) 


0.541(15) 


0.797(18) 


6 


1.251 


1.201(66) 


0.554(21) 


0.771(33) 


0.763(22) 


0.540(17) 


0.797(19) 


7 


1.121 


1.202(71) 


0.551(22) 


0.770(35) 


0.763(24) 


0.539(18) 


0.796(20) 


8 


1.017 


1.202(76) 


0.548(24) 


0.770(37) 


0.762(26) 


0.538(19) 


0.795(22) 



Table 3: Perturbative matching (cf. Eq. ()2.29p ) for various choices of the matching scale 



whose form is motivated by the perturbative series, with coefficients 

si = 7oln2, (3.3) 
1 



S2 = 71 In 2 + 



^(70)^ + ^070 



(In 2)2 . (3.4) 



It is worth stressing that si is universal and independent of A'^f, whereas S2 carries a 
dependence upon A^f via 60 and 71, with the latter coefficient introducing a scheme 
dependence. In our fits we truncated the polynomial at 0{u^). The fits have been 
performed with si fixed to its perturbative value and S2, S3 left as free parameters. 
Fit results are shown in Figs. [3l [5] and El Fitted values of S2 turned out to be close 
to the perturbative prediction of Eq. (j3.4p . with the exception of Q'^. 

Once the continuous SSFs have been obtained as functions of the renormalised 
coupling, the ratios c (cf. Eq. (|2.16|) ) are obtained recursively. The low-energy scale 
A^had = -^max is implicitly defined in this work through the condition gsp(Linax) = 
4.61, as explained in [18]. This scale is chosen so that the renormalisation con- 
stants Z((7o, a^Uhad) can be computed in the accessible ^o-range commonly used in 
large-volume simulations. The non-perturbative RG running of the six operators of 
interest are shown in Fig. [SI 

As discussed in our former quenched study [4], the main criterion for selecting 
robust schemes amounts to checking that the systematic uncertainty present in our 
final results, due to the NLO truncation of the perturbative matching at the scale 
Upt = 2"'//had) is well under control. This in turn requires an estimate of the size of 
the NNLO contribution to c. To this purpose we have re-computed c with two dif- 
ferent ansatze for the NNLO anomalous dimensions 72: (i) we set I72/71I = I71/70I) 
(ii) we perform a two-parameter fit to the SSF with si,S2 fixed to their perturba- 
tive values and 53,54 left as free parameters, and then estimate 72 by equating the 
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resulting value of S3 to its perturbative expression 



S3 = 72 In 2 + [7071 + 26o7i + 6170] (In 2)^+ 



+ g7o + &o7o + 3&070 (In 2)3 . 



(3.5) 



The optimal schemes specified in Tabled] are precisely those for which the aforemen- 
tioned determinations of the effective 72 lead to the smallest discrepancies between 
the corresponding universal factors c. 

The effect of varying the perturbative matching point in the optimal schemes 
is described by Table [3l We see that numbers are very stable for n > 6, while the 
uncertainty increases with n due to progressive error accumulation at each step. 
Final results, reported in the second column of Table [7| refer to n = 7. Note that 
typical relative errors are as big as 5%, which may result in a sizeable error in 
hadronic matrix elements, solely due to renormalisation. 

4 Connection to hadronic observables 

Having computed the universal evolution factors c(//had)) which provide the RG- 
running from the low energy matching scale ^had to a formally infinite one, we 
proceed to establish the connection between bare lattice operators and their RGI 
counterparts. The latter, defined in Eq. ()2.16p from the integration of the Callan- 
Symanzik equation, are related to the bare operators used in lattice simulations via 
a total renormalisation factor Z^Qi{go), defined as 



The Zrgi factor does not depend on any renormalisation scale and carries a depen- 
dence upon the renormalisation condition only via cutoff effects. 

In order to obtain Z(go, O/Uhad)) we follow [18] and compute Z{go,afi) at three 
values of the lattice spacing, namely /5 = {5.20, 5.29, 5.40}, which belong to a range 
of inverse couplings commonly used for simulations of two-flavour QCD in physically 
large volumes. Simulation parameters and results are collected in Table H] for the 
relativistic operators Qf and in Table [5] for the static ones Q'^ ■ 

While the simulation at (/3 = 5.20, L/a = 6) is exactly at the target value 
for 5gF(Lmax), corresponding to Z{go,a^ha.d), the simulations at the other /3 values 
require a slight interpolation. We adopt a fit ansatz, motivated by Eq. (j2.16p . 



in order to interpolate the Z factors between the values of g^p straddling the target 
value 5sp(-^^max) = 4.61. Note that the fits take into account the (independent) 



Zrgi{9o) = ^(90,a/^had)c(/Uhad) • 



(4.1) 



ln(Z) = ci +C2ln(g|p) , 



(4.2) 
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Z+ 


z- 


5.20 


0.13600 


4 


3.65(3) 


0.7547(19) 


0.4797(12) 






6 


4.61(4) 


0.7715(20) 


0.4383(11) 


o.zy 


n 1 1 

U. lOD'il 






U. ( OOol^i ( ) 


U.OU t\J\LL) 






6 


4.297(37) 


0.7749(24) 


0.4644(13) 






8 


5.65(9) 


0.8036(26) 


0.4339(12) 


5.40 


0.13669 


4 


3.188(24) 


0.7591(16) 


0.5342(11) 






6 


3.864(34) 


0.7709(21) 


0.4871(13) 






8 


4.747(63) 


0.7938(22) 


0.4583(11) 



Table 4: Results for Z+ and Z with ct set to its 2-loop value. The values of g^^ are 
from [26]. The hopping parameters k, used in the simulations are the critical ones (wcr) 
of [27]. 



/3 




L/a 




zw 


2(2) 


2(3) 


2(4) 


5.20 


0.13600 


4 
6 


3.65(3) 
4.61(4) 


0.7793(17) 
0.7118(17) 


0.9741(16) 
0.9409(15) 


0.8681(16) 
0.7857(13) 


0.8317(13) 
0.7921(13) 


5.29 


0.13641 


4 
6 
8 


3.394(17) 
4.297(37) 
5.65(9) 


0.7862(16) 
0.7275(20) 
0.6612(19) 


0.9766(16) 
0.9431(18) 
0.9150(16) 


0.8768(15) 
0.7992(16) 
0.7337(15) 


0.8397(14) 
0.8017(15) 
0.7619(14) 


5.40 


0.13669 


4 
6 
8 


3.188(24) 
3.864(34) 
4.747(63) 


0.7972(15) 
0.7378(18) 
0.6840(16) 


0.9805(14) 
0.9434(17) 
0.9231(14) 


0.8864(14) 
0.8098(15) 
0.7529(13) 


0.8497(12) 
0.8094(14) 
0.7781(12) 



Table 5: Results for Z'^^^ with ct set to its 2-loop value. The values of gl-^ are from [26]. 
The hopping parameters k used in the simulations are the critical ones (kci) of [27]. 



errors of both Z and g"^^. Moreover, we have conservatively augmented the fit errors 
by the difference between the fit results of Eq. (14. 2|) and the results from a naive 
two-point linear interpolation in ^|p. The coefficients C2 of the fits (|4.2|) deviate in a 
range of 7% — 30% from the lowest order coefficients 7o/(26o), signalling the presence 
of moderate higher-order perturbative effects. 

The resulting numbers for the renormalisation factors at the low energy match- 
ing scale, and also for the RGI renormalisation factors .^rgi(5o)) are collected in 
Table m The first error of the stems from the error of Z factors, whereas the 

second accounts for the uncertainties in the universal factors c. Note that only the 
first of these errors should be added in quadrature to the error of the bare hadronic 
matrix elements, once these become available from future computations, in order 
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•^RGI 


5.20 
5.29 
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0.7715(20) 
0.7825(27) 
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0.927(2)(55) 
0.940(3) (56) 
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0.4383(11) 
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0.4623(25) 


0.241(1)(10) 
0.251(1)(10) 
0.255(1)(10) 




zw 


7(1) 

^RGI 


^(2) 


7(2) 

^RGI 


5.20 
5.29 
5.40 


0.7118(17) 
0.7093(27) 
0.6904(40) 


0.548(1)(28) 
0.546(2) (28) 
0.532(3)(28) 


0.9409(15) 
0.9374(30) 
0.9233(46) 


0.718(1)(26) 
0.715(2)(26) 
0.704(4) (26) 


/3 


^(3) 


^(3) 


zw 


^RGI 


5.20 
5.29 
5.40 


0.7857(13) 
0.7836(40) 
0.7567(75) 


0.423(1)(15) 
0.422(2)(15) 
0.408(4) (14) 


0.7921(13) 
0.7916(18) 
0.7807(38) 


0.631(1)(17) 
0.630(1)(17) 
0.621(3)(17) 



Table 6: Results for Z+, Z , Z^^^ and Z^^-^, ^rgh -Z^rgi for three bare gauge coupling 
values corresponding to our low-energy matching point at gl^ — 4.61 in the SF scheme. 

to obtain the total error of the renormalised quantity, at a given lattice spacing. 
The second error, which is entirely unrelated to the discretisation of the theory, 
should only be added in quadrature to the continuum extrapolated hadronic matrix 
element. For the sake of convenience, a representation of the numerical results of 
Table El by interpolating polynomials is also adopted, i.e. 

Zkgi = ao + ai{f3- 5.2) + a2(/3 - 5.2)2 ^ (4 3) 

which can be used at any intermediate value of /5 between /3 = 5.20 and /? = 5.40. 
Fit coefficients are reported in Table [7] for the various operators. The uncertainty of 
the RGI constants at intermediate points may be easily obtained from those at the 
simulation points, see Table [U by linear interpolation. 

As a final remark, we observe that the simulation of the renormalisation factors 
at /3 = 5.20, L/a = 4 is not at the target value for g^p. We used is as a check of 
the independence of the Z^qi , computed via Eq. (|4.1|) from the low energy matching 
scale. Specifically, the two measured values of Z-factors at /? = 5.20 have been 
used in order to extrapolate the renormalisation constants at 5sF(-^^max/2) = 3.0318, 
where the non-perturbative matching with the universal evolution factors c has been 
subsequently performed. Results turned out to be fully compatible with those quoted 
in Table El 
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n 1 

(XI 


0,2 


Qt 


1.202(71) 


0.9270 


0.1741 


-0.2973 


Qi 


0.551(22) 


0.2414 


0.1431 


-0.3853 


Q't 


0.770(35) 


0.5481 


0.0285 


-0.5546 


Q't 


0.763(24) 


0.7179 


0.0010 


-0.3407 


Q't 


0.539(18) 


0.4235 


0.0411 


-0.5962 


Q't 


0.796(20) 


0.6305 


0.0291 


-0.3723 



Table 7: Universal factors c and coefficients of the interpolating polynomials of the RGl 
renormalisation constants, see Eq. ()4.3p . Uncertainties are discussed in the text. 

5 Conclusions 

Using standard SF methods, we have performed a fully non-perturbative compu- 
tation of the renormalisation and RG running of several four-fermion operators in 
= 2 QCD. We have considered the two operators made of four relativistic quark 
fields with a left-left Dirac structure and the complete basis of operators with two 
static and two relativistic quarks. The Wilson lattice actions have been implemented 
for both the gauge and the fermionic parts, the latter with a non-perturbatively 
tuned Clover term. The HYP2 discretisation of the static quark turned out to be 
the less noisy choice, after comparison with other options. Only the parity-odd parts 
of the operators have been analysed, as their renormalisation pattern is unaffected 
due to the loss of chiral symmetry by the regularisation. 

Our results are an essential building block for any A'^f = 2 computation of 
quantities like Bk and Bb- Nevertheless, their precision is somewhat limited by 
increased statistical fluctuations at the three strongest couplings and by the lack of 
a fourth, finer, lattice resolution which would improve the continuum extrapolation 
of the operator SSFs. This could lead to a potentially unsatisfactory total error 
on hadronic matrix elements. Future refinement (besides using a two-loop estimate 
of Ct throughout the runs and increased statistics at the three strongest couplings) 
is necessary, either by simulating closer to the continuum limit, or by completely 
removing leading order discretisation effects from the simulations. 
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Figure 2: Continuum extrapolation of the SSFs for (left) and (right). The renor- 
malised coupling increases from top to bottom. Blue dotted lines and the blue cross at 
a/L = correspond to weighted averages of the L/a = 8, 12 data, red dashed lines and the 
red a/L = open point to linear extrapolations in a/L of the three data. 
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Figure 3: The step-scaling functions (t+ and (discrete points) as obtained non- 
perturbatively. The shaded area is the one sigma band obtained by fitting the points to 
a polynomial as discussed in the text. The dotted (dashed) line is the LO (NLO) perturba- 
tive result. 
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Figure 4: Continuum extrapolation of the SSFs for Q'^ (left) and Q'J (right). The renor- 
malised coupling increases from top to bottom. Blue dotted lines and the blue cross at 
a/L = correspond to weighted averages of the L/a = 8, 12 data, red dashed lines and the 
red a/L — open point to linear extrapolations in a/L of the three data. 




Figure 5: The step-scaling functions cr^^^ and cr'-^-' (discrete points) as obtained non- 
perturbatively. The shaded area is the one sigma band obtained by fitting the points to 
a polynomial as discussed in the text. The dotted (dashed) line is the LO (NLO) perturba- 
tive result. 
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Figure 6: Continuum extrapolation of the SSFs for Q'J (left) and Q'4 (right). The renor- 
malised coupling increases from top to bottom. Blue dotted lines and the blue cross at 
a/L = correspond to weighted averages of the L/a = 8, 12 data, red dashed lines and the 
red a/L — open point to linear extrapolations in a/L of the three data. 
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Figure 7: The step-scaling functions cr^'^^ and cr'-'*-' (discrete points) as obtained non- 
perturbatively. The shaded area is the one sigma band obtained by fitting the points to 
a polynomial as discussed in the text. The dotted (dashed) line is the LO (NLO) perturba- 
tive result. 
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Figure 8: RG-running of of the four-quark operators obtained non-perturbatively (discrete 
points) at specific values of the renormalisation scale fi, in units of A. The lines are per- 
turbative results at the order shown for the Callan-Symanzik /3 function and the operator 
anomalous dimension 7. 
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Table 8: Numerical values of the renormalisation constants Z+, Z~ and the step scaling 
functions E"'", E~ at various rcnormaliscd SF couplings and lattice spacings. Data at g^p = 
0.9793, 1.1814, 1.5031 have been obtained with ct evaluated in one-loop perturbation theory. 
The remaining data have been obtained with ct evaluated in two-loop perturbation theory. 
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Table 9: Numerical values of the renormalisation constants 2^^\ Z^^^ and the step-scaling 
functions S^'^-', S^^^ with HYP2 action at various renormalised SF couplings and lattice 
spacings. Data at = 0.9793, 1.1814, 1.5031 have been obtained with ct evaluated in 
one-loop perturbation theory. The remaining data have been obtained with ct evaluated in 
two-loop perturbation theory. 
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Table 10: Numerical values of the renormalisation constants Z^^\ Z^^^ and the step-scaling 
functions S^^\ S^''^ with HYP2 action at various renormalised SF couplings and lattice 
spacings. Data at g^p = 0.9793, 1.1814, 1.5031 have been obtained with ct evaluated in 
one-loop perturbation theory. The remaining data have been obtained with ct approximated 
in two-loop perturbation theory. 
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